Warm up example on Irreps

Consider the tensor product of two random single qubit unitaries

Let' s print a random single qubit unitary :        i.e. a random U in U(2)

Urand = RandomVariate[CircularUnitaryMatrixDistribution[2]];
Urand // MatrixForm

Let' s print the tensor product of two random single qubit unitaries:       i.e.  U  \otimes U  for a random  U in U(2)

KroneckerProduct[Urand, Urand]//MatrixForm

This looks a mess right ...  But what if we rotate into the Bell basis?   (Remember this is done by applying a Hadamard and a CNOT) 

Let' s define these gates :

CNOT = {{1,0,0,0},{0,1,0,0},{0,0,0,1},{0,0,1,0}};
CNOT // MatrixForm

Hadamard = (1/(2^(1/2))){{1,1},{1,-1}};
Hadamard // MatrixForm

Now let' s rotate the random unitary into this basis - we see that it is now BLOCK DIAGONAL

(KroneckerProduct[Hadamard,IdentityMatrix[2]]).CNOT.KroneckerProduct[Urand, Urand].CNOT.(KroneckerProduct[Hadamard,IdentityMatrix[2]])  // MatrixForm

This wasn't a freak for that example - it works for any random single qubit unitary

Urand = RandomVariate[CircularUnitaryMatrixDistribution[2]];
KroneckerProduct[Urand, Urand]//MatrixForm
(KroneckerProduct[Hadamard,IdentityMatrix[2]]).CNOT.KroneckerProduct[Urand, Urand].CNOT.(KroneckerProduct[Hadamard,IdentityMatrix[2]])  // MatrixForm

That is, we have taken a tensor product representation and reduced it to a direct sum representation. Equivalently, we have block diagonalized the representation.

We can play the same game for a random unitary in SU(2) instead of U(2)

Don’t worry too much about the details - I'm doing the same as above but just making sure that the random unitary has determinant 1

Urand = RandomVariate[CircularUnitaryMatrixDistribution[2]];
Urand //MatrixForm;
Correction = (Conjugate[Det[Urand]] IdentityMatrix[2])^(1/2);
SUrand = Correction.Urand;
SUrand // MatrixForm
Round[KroneckerProduct[SUrand , SUrand ],0.001]//MatrixForm
BellBasisOutput = (KroneckerProduct[Hadamard,IdentityMatrix[2]]).CNOT.KroneckerProduct[SUrand ,SUrand].CNOT.(KroneckerProduct[Hadamard,IdentityMatrix[2]]);
Round[BellBasisOutput,0.001]  // MatrixForm

Twirling example (Ignore this for now)

del={{4}};
SymOutput=Delete[BellBasisOutput,del];
SymOutput=Delete[Transpose[SymOutput],del];
SymOutput=Transpose[SymOutput];
Round[SymOutput, 0.001] // MatrixForm
Round[Det[SymOutput], 0.001]

randoutput[X_] :=Module[ {Urand},Urand= RandomVariate[CircularUnitaryMatrixDistribution[2]];KroneckerProduct[Urand, Urand].X.ConjugateTranspose[KroneckerProduct[Urand, Urand]]]

rho= {{0.2,0,0,0},{0,0.3,0,0},{0,0,0.1,0},{0,0,0,0.4}};

rho = RandomVariate[CircularUnitaryMatrixDistribution[4]].rho.ConjugateTranspose[RandomVariate[CircularUnitaryMatrixDistribution[4]]];
rho//MatrixForm

averagedstate = Round[Mean[Table[randoutput[rho], {i, 1, 1000000} ]],0.001];
averagedstate // MatrixForm

RotateToBell[rho_]:= (KroneckerProduct[Hadamard,IdentityMatrix[2]]).CNOT.rho.CNOT.(KroneckerProduct[Hadamard,IdentityMatrix[2]])

RotateToBell[averagedstate]// MatrixForm

SymProj = {{1,0,0,0}, {0,(1/2),(1/2),0}, {0,(1/2),(1/2),0}, {0,0,0,1}};

AntiProj =  {{0,0,0,0}, {0,(1/2),((-1)/2),0}, {0,((-1)/2),(1/2),0}, {0,0,0,0}};

Tr[rho.SymProj]

Tr[rho.AntiProj]

IdSym = {{1,0,0,0}, {0,1,0,0}, {0,0,1,0},{0,0,0,0}};

IdAnti =  {{0,0,0,0}, {0,0,0,0}, {0,0,0,0},{0,0,0,1}};

Round[(1/3)Tr[rho.SymProj]IdSym+ Tr[rho.AntiProj]IdAnti , 0.001]// MatrixForm